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A  methodology  is  presented  to  sample  efficiently  configurations  of  reacting  mixtures  in 
the  reaction  ensemble  Monte  Carlo  simulation  technique.  A  cavity-biasing  scheme  is  used, 
which  more  effectively  samples  configurations  than  conventional  random  sampling.  Akin 
to  other  biasing  schemes  that  are  implemented  into  insertion-based  Monte  Carlo  methods 
such  as  Gibbs  ensemble  Monte  Carlo,  the  method  presented  here  searches  for  space  in 
the  reacting  mixture  whereby  the  insertion  of  a  product  molecule  is  energetically  favoured. 
This  sampling  bias  is  then  corrected  in  the  acceptance  criteria.  The  approach  allows  for  the 
study  of  reacting  mixtures  at  high  density  as  well  as  for  polyatomic  molecular  species.  For 
some  cases,  the  method  is  shown  to  increase  the  efficiency  of  the  reaction  steps  by  a  factor 
greater  than  20.  The  approach  is  shown  to  be  readily  generalized  to  other  biasing  schemes 
such  as  orientational-biasing  of  polar  molecules  and  configurational-biasing  of  chain-like 
molecules. 


1.  Introduction 

The  reaction  ensemble  Monte  Carlo  simulation  method 
(RxMC)  is  a  powerful  simulation  tool  for  studying 
the  behaviour  of  chemical  reaction  equilibria  [1,  2]. 
The  RxMC  technique  can  be  used  to  predict  the  shift 
in  reaction  equilibria  for  an  ideal-gas  phase  reaction  in 
a  non-ideal  environment  such  as  a  condensed  phase, 
a  solvent,  or  a  porous  material.  Recent  applications  of 
the  RxMC  method  include  reactions  of  plasmas  [3], 
reactions  in  porous  carbons  [4-7]  and  porous  mem¬ 
branes  [8],  reactions  under  shock  conditions  [9],  and 
reactions  in  supercritical  fluid  solvents  [10,  11]. 

Comprehensive  reviews  of  the  RxMC  method  [12]  and 
its  applications  to  date  [13]  can  be  found  elsewhere. 

As  part  of  the  configuration  sampling  algorithm  in  the 
RxMC  method,  the  insertion  and  deletion  of  molecules 
is  required  to  satisfy  the  reaction  equilibrium  condition 
Vjiiij  —  0,  where  v/7  is  the  stoichiometric  coefficient 
of  species  i  in  chemical  reaction  j,  /x,-  is  the  chemical 
potential  of  species  i  in  the  reacting  mixture,  and  C  is 
the  total  number  of  reactants  and  products.  Analogous 
to  other  Monte  Carlo  (MC)  methods  such  as  grand 
canonical  ensemble  MC  (GCMC)  [14,  15]  and  Gibbs 
ensemble  MC  (GEMC)  [16,  17]  that  require  the  insertion 
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of  molecules  into  the  simulation  box,  random  insertion 
fails  for  high  density  systems  or  for  multi-site  models 
such  as  polymers  in  even  moderately  dense  systems. 
Conventional  random  sampling  fails  due  to  the  high 
probability  of  the  inserted  molecule  encountering 
an  overlap  with  molecules  already  in  the  simulation 
box,  generating  a  configuration  that  is  energetically 
unfavorable  and  most  probably  rejected. 

Previous  applications  of  the  RxMC  method  success¬ 
fully  used  conventional  random  sampling  because  low 
density  phases  were  studied  or  simple  molecular  models 
were  considered.  In  this  work,  we  present  a  cavity-bias 
sampling  method  for  the  RxMC  technique  that  makes 
it  possible  to  simulate  reaction  equilibria  at  conditions 
where  conventional  sampling  fails.  The  method,  termed 
cavity-bias  reaction  ensemble  Monte  Carlo  (CB-RxMC) 
is  analogous  to  bias  sampling  schemes  used  in  the 
GCMC  [18,  19]  and  GEMC  simulation  methods  [20], 
As  the  name  implies,  the  method  searches  for  cavities  in 
the  system  into  which  product  molecules  can  be  success¬ 
fully  inserted.  The  method  can  be  easily  generalized  to 
include  other  biasing  methods  such  as  orientational- 
biasing  and  configurational-biasing. 

In  the  following,  we  briefly  review  the  essentials  of 
the  RxMC  method,  after  which  we  present  the  formal¬ 
ism  of  our  cavity-bias  sampling  method  for  the  RxMC 
technique.  We  present  an  example  of  the  method  to 
demonstrate  its  efficiency. 
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2.  Formalism 

2.1.  Reaction  ensemble  Monte  Carlo 

The  RxMC  method  is  designed  to  minimize  the  Gibbs  free 
energy  thus  determining  the  true  chemical  equilibrium 
state  irrespective  of  rate  limitations.  RxMC  requires 
intermolecular  potentials  for  the  molecular  species  that 
are  present  in  the  reacting  mixture.  RxMC  also  requires 
inputting  the  ideal-gas  internal  modes  for  each  react¬ 
ing  species,  information  which  is  readily  available 
in  standard  sources  [21,  22]  or  can  be  generated  using 
quantum  mechanics  methods.  Finally,  the  particular 
reactions  occurring  in  the  system  must  be  specified. 

Implementation  of  RxMC  provides  information 
about  the  chemical  equilibrium  state,  such  as  the  density 
of  the  reacting  mixture,  mole  fractions  of  the  reacting 
species,  the  change  in  the  total  number  of  moles,  and 
the  internal  energy.  RxMC  can  be  performed  in  many 
different  types  of  ensembles  [12],  including  canonical, 
isothermal-isobaric,  Gibbs  and  other  less  common 
ensembles  [23].  Furthermore,  RxMC  can  be  performed 
for  multiple  reactions  and  multiple  phases  [24].  RxMC 
does  not  simulate  bond  breaking  or  formation;  rather, 
RxMC  directly  samples  forward  and  reverse  reaction 
steps  as  Monte  Carlo-type  moves  according  to  the 
stoichiometry  of  the  reactions  being  considered. 

The  isothermal-isobaric  ensemble  version  of  the 
RxMC  method  for  a  set  of  linearly  independent  chemi¬ 
cal  reactions  involves  the  following  trial  moves: 

(1)  a  change  in  the  position  or  orientation  of  a  molecule 
which  is  chosen  at  random; 

(2)  a  random  change  in  the  simulation  box  volume;  and 

(3)  for  a  randomly  chosen  reaction  and  randomly 
chosen  reaction  direction  (forward  or  reverse), 
reactant  molecules  are  randomly  chosen  and  deleted 
while  product  molecules  either  replace  reactant 
molecules  or  are  inserted  randomly  into  the 
simulation  box. 

Step  (1)  ensures  that  thermal  equilibrium  is  established 
for  the  specified  temperature,  Step  (2)  ensures  that 
mechanical  equilibrium  is  established  for  the  specified 
pressure,  and  Step  (3)  ensures  that  chemical  equilibrium 
is  established  for  the  specified  reactions. 

Step  (3)  requires  the  insertion  of  product  molecules 
into  the  simulation  box.  For  reactions  where  the  number 
of  moles  decreases  or  remains  unchanged,  the  product 
molecules  can  simply  be  placed  in  the  space  previously 
occupied  by  the  deleted  reactant  molecules.  For  reac¬ 
tions  where  the  number  of  moles  increases,  some  of  the 
product  molecules  can  replace  reactant  molecules  but 
the  remaining  product  molecules  must  be  inserted 
into  the  simulation  box.  In  the  conventional  sampling 


scheme  this  entails  randomly  choosing  a  position  in  the 
box.  For  medium  to  high  density  systems,  this  step  can 
result  in  a  prohibitively  low  acceptance  rate  since  a 
random  insertion  will  almost  always  result  in  an  overlap 
with  molecules  already  in  the  box.  Such  an  overlap  is 
energetically  unfavourable  and  therefore  the  attempted 
reaction  step  will  nearly  always  be  rejected.  Below  we 
describe  an  algorithm  that  allows  us  to  bias  the  insertion 
of  product  molecules  in  such  a  way  that  cavities  in  the 
system  are  found.  This  method  of  insertion  would  bias 
the  simulations  if  the  ordinary  acceptance  rules  were 
used.  Therefore,  we  next  demonstrate  that  the  correct 
distribution  of  configurations  can  be  sampled  if  the 
acceptance  rule  for  this  step  is  modified  appropriately. 


2.2.  General  approach  for  bias 
Monte  Carlo  sampling 

For  completeness,  we  first  consider  the  general  approach 
for  bias  Monte  Carlo  sampling  (see  [25]  for  a  more 
thorough  account).  Let  K(o  — »■  n)  be  the  flow  of  con¬ 
figuration  o  — >  n,  where 


K(o  — >  n)  —  N(o)  x  a(o  -»  «)  x  acc{o  -*  n).  (1) 


N(o )  is  the  probability  of  being  in  configuration  o ; 
a(o  -»  n)  is  the  probability  of  generating  configuration  n; 
and  acc(o  — >  n)  is  the  probability  of  accepting  the  move 
from  o  ->  n.  N,  often  referred  to  as  the  probability 
density  or  distribution  of  configurations,  depends  on 
the  details  of  the  ensemble  and  can  be  determined  from 
the  partition  function.  For  example,  the  probability  of 
finding  configuration  rN  in  the  canonical  ensemble  is 


N(  rN)  = 


exp[-/ff/(rN)] 

/  drNexp  [— /lC/(rN)] 


(2) 


where  t/(rN)  is  the  configurational  energy  of  N  particles 
and  f—  l/kB?"  (kB  is  the  Boltzmann  constant  and  T  is 
the  temperature). 

To  guarantee  a  correct  sampling  scheme,  the  detailed 
balance  condition  is  imposed 


K(o  -*  n)  —  K(n  — ►  o) 


(3) 


or  more  specifically 


N(o)  x  a(o  -*  n)  x  acc(o  — »■  n) 

—  N(n)  x  a(n  — >  o)  x  acc(n  —>  o).  (4) 
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In  the  original  Metropolis  scheme,  a  is  chosen  to  be  a 
symmetric  matrix,  i.e.  a(o  — >  n)  —  a(n  — >  o).  For  a 
biased  Monte  Carlo  algorithm,  we  wish  to  generate 
configurations  in  a  biasing  manner  thus  making 
a(o  ->  n)  /  a(n  — »■  o).  It  is  clear  from  equation  (4), 
that  since  a(o  — »■  n)  f  a(n  ->  o)  and  7V(o)  and  7V(«)  are 
predetermined  we  must  correct  for  biasing  a  by 
changing  the  ratio  accip  — >  n)/acc{n  —>  o). 

Suppose  that  we  have  developed  a  Monte  Carlo 
scheme  to  generate  trial  configurations  with  a  prob¬ 
ability  that  depends  on  the  potential  energy  of  that 
configuration 

a(o  n)  =  / [77(h)]  (5) 


qmli  is  the  quantum  partition  function  for  the  internal 
modes  of  an  isolated  molecule  of  species  i,  which 
includes  vibrational,  rotational,  and  electronic;  A,  is  the 
thermal  de  Broglie  wavelength  of  species  i;  and  V  is 
the  total  volume  of  the  system  [1,  2].  Equation  (9)  is 
appropriate  for  both  forward  and  reverse  reaction  steps 
(§y  =  1  for  a  forward  step  and  £,  =  —  1  for  a  reverse  step), 
where  the  stoichiometric  coefficients  are  taken  to  be 
positive  for  product  species  and  negative  for  reactant 
species.  U,  is  the  configurational  energy  between 
molecule  i  and  all  other  molecules  in  the  simulation 
box,  so  for  a  deleted  reactant  molecule  Ufn)  —  0,  while 
for  an  inserted  product  molecule  I7,(o)  =  0.  From 
equation  (9)  we  can  write  the  acceptance  rule  for  a 
reaction  step  using  conventional  sampling  in  RxMC  as 


and  while  for  the  reverse  move 

ot{n  -*  o)  —f[U{o)\. 
From  equation  (4)  then 


(6) 


acc(o  — »■  n)  —  min 


i-n 

7=1 


Nil 

(N,+  Vj£j)l\  a; 


x  exp [-P(Uj(n)  -  77, (o))] 


(10) 


acc(o  n)  _  f[U{oj]  N(ri) 

acc(n  -*  o)  ~  /'  [[/(/?)]  Nip) '  7  Similar  to  other  biasing  algorithms  [25],  we  can  define 

a  Rosenbluth  factor 


A  possible  acceptance  rule  that  obeys  this  condition  is 


acc{o 


n)  —  min 


/  f[U(o)]N(n)\ 
V  \f[U(n)]N(o))- 


(8) 


k, 

wi(n )  =  exP l-PUi(l )], 


/=  i 


(11) 


It  is  clear  then  that  we  can  introduce  an  arbitrary 
biasing  function /[  77]  in  the  sampling  scheme,  provided 
that  the  acceptance  rule  is  modified  in  such  a  way  that 
the  bias  is  removed  [25].  Note  that  the  same  biasing 
function  must  be  used  to  generate  the  o  configurations 
and  the  n  configurations. 


which,  for  a  single  attempted  move,  is  determined  by 
generating  k,  trial  configurations  for  species  i.  i.e.  when  i 
is  a  reactant  molecule  kt  trial  deletions  and  when  i  is  a 
product  molecule  kj  trial  insertions.  From  these  trial 
configurations  we  can  select  a  configuration,  m,  for 
species  i  with  a  probability 


2.3.  Derivation  of  acceptance  rule  for  CB-RxMC 

For  a  reaction  step  in  RxMC,  the  ratio  of  the 
probability  densities  of  the  old  ‘o’  and  the  new  V 
configurations  can  be  written 

Njn)  =  A  Ay!  /qinuiV 
N(o)  )J(M  +  C,?/)!  V  A?  ) 

x  exp [-P(Ui(ri)  -  Ui(o))\  (9) 


(  =  exp[-^7/,(m)] 

Etiexp[-i8^(/)] 


(12) 


where  f(U,)  is  a  biasing  function  that  ensures  energeti¬ 
cally  favourable  configurations  are  generated.  The 
expression  for  the  complete  reaction  step  is  then 


9 


9 


f\U{o)\  =  Y\f[U,(o)\  =  l 

7=1  7=1 


exp[-;6 77,(0)] 

Et,  exP[-/17 /?(/)] 


(13) 


where  Cj  is  the  total  number  of  species  in  reaction  /; 
Nj  is  the  total  number  of  molecules  of  species  i;  vyV  is 
the  stoichiometric  coefficient  of  species  i  in  reaction  /; 
is  the  molecular  extent  of  reaction  for  reaction  /; 


where  we  can  write  an  analogous  expression  for  f[U(n)\. 

Substituting  expressions  for  /[77(o)],  ,/[7/(«)]  and 
equation  (9)  into  equation  (7)  we  get  the  probability 
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of  accepting  the  reaction  step  attempt  as 


acc(o  -»  n ) 
acc(n  — >■  o) 


n 


W! 

^int,  i  ^  \ 

(Nj  +  Vj£j)\ 

l  A?  J 

exp[— /3(t/,(«)  -  Uj(o))\ 


exp[-j RUj(o)\ 

Eti  exp [~pUf(l)]  exp [~pUi(n)] 


(14) 


while  after  simplifying  we  get 

acc(o  —>  n) 
acc(n  —>  o) 


r  n^ 

^int,  i  V  ^ 

exp[— /i !/*(/)]" 

(Nj  +  Vj£j)\ 

l  A  ?  ) 

Eti  exp [~pUf(l)\_ 

(15) 

Finally  then,  the  acceptance  rule  for  a  reaction  step 
that  satisfies  the  detailed  balance  condition  is 


acc(o  — »■  n)  —  min 


'•n 

7=1 


Nil 


( V,  +  Vj&)l 


E/Ii  exp[— /!£/”(/ )] 
Eti  exp[-/JC/,°(/)] 


(16) 


Recall  that  when  i  is  a  deleted  reactant  molecule, 
U"{1)  —  0  for  all  /’ s,  so  Etii  exp[— ySC/f(/)]  =  kj. 

Likewise,  when  i  is  an  inserted  product  molecule, 
£/?(/)  =  0  for  all  /’ s,  so  E/=i  exp[— /)£/?(/ )]  =  kj. 

Finally,  note  that  when  k,  =  1  for  all  species, 
equation  (16)  reduces  to  the  conventional  sampling 
acceptance  rule  given  in  equation  (10). 

With  minor  extensions,  the  above  scheme  can  be 
applied  to  include  orientational-biasing  algorithms  for 
molecular  models  with  a  strong  dependence  on  relative 
orientation  (e.g.  polar  and  hydrogen-bonded  molecules 
and  multi-atomic  models)  or  configurational-biasing 
schemes  for  chain-like  molecules  [25],  Using  the 
Rosenbluth  factors,  equation  (15)  can  be  generalized 
to  accommodate  such  biasing  schemes  by  writing 


acc(o  —>  n)  y-f 

m 

^Wi(n) 

acc(n  —*■  o)  || 

(N  +  r/, $/)! 

{  A'  J 

Wj(o) 

(17) 


Appropriate  evaluation  of  the  Rosenbluth  factors  for 
the  particular  scheme  is  required. 

2.4.  CB-RxMC  algorithm 

Next,  we  provide  a  detailed  outline  for  implementing  the 
CB-RxMC  algorithm.  For  simplicity,  consider  a  system 
of  J  reactions  where  all  reacting  species  are  modelled 
as  single  spheres  and  electrostatic  contributions  are 
ignored.  (Again,  this  algorithm  can  be  readily  extended 
to  include  multi-atomic  and  electrostatic  molecular 
models.) 

Step  1 :  Randomly  choose  reaction  j  and  the  reaction 
direction  (forward  or  reverse). 

Step  2:  Randomly  select  reactant  molecules  and 
determine  exp[— /)[/?( 1)]  for  each  of  these 
reactant  molecules,  where  /  =  1 . 

Step  3:  For  each  reactant  molecule,  generate  k —  1 
trial  configurations  denoted  bi  ■  ■  ■  bk,  and 
determine  E/=2  exp[— PU°(b/)\.  Note  that  for 
reactants  whose  centre-of-mass  will  be 
replaced  by  the  centre-of-mass  of  a  product 
molecule,  k,  —  1 . 

Step  4:  For  product  molecules  whose  centre-of-mass 
will  replace  the  centre-of-mass  of  a  deleted 
reactant  molecule,  determine  exp[— PU"(  1)] 
where  /  =  1,  and  set  kj  —  1. 

Step  5:  For  the  remaining  product  molecules  that 
need  to  be  inserted  for  reaction  j,  generate  kj 
trial  configurations,  denoted  b\  . . .  bk„  by 
randomly  inserting  the  product  molecule 
into  the  simulation  box  and  determine 
Efii  exp[— pU"(bi)].  From  these  kt  insertions, 
select  one  with  the  probability  given  in 
equation  (12). 

Step  6:  The  attempted  reaction  step  for  reaction  j 
is  accepted  with  the  probability  given  in 
equation  (16). 

It  is  implicitly  assumed  that  when  product  molecules 
are  inserted  they  are  placed  in  the  locations  previously 
occupied  by  reactant  molecules.  This  is  not  a  require¬ 
ment  of  the  RxMC  method  but  for  most  cases  will 
enhance  the  phase  sampling  efficiency.  Moreover, 
replacing  reactant  molecules  with  product  molecules 
which  have  similar  molecular  model  character  (e.g.  size, 
atomistic  detail,  or  polarity)  will  further  increase 
sampling  efficiency.  However,  the  reactant-to-product 
replacement  must  be  consistent  for  the  forward  and 
reverse  directions  for  a  particular  reaction  throughout 
the  simulation  run.  For  example,  if  a  type  a  molecule  is 
replaced  by  a  type  b  molecule  in  the  forward  direction  of 
a  reaction  step,  then  a  type  b  molecule  must  be  replaced 
by  a  type  a  molecule  in  the  reverse  direction  [1,  2], 
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An  alternative  but  less  efficient  approach  (for  most 
cases)  would  be  to  randomly  insert  all  product 
molecules.  The  necessary  minor  adjustments  to  the 
above  algorithm  would  be  to  ignore  the  last  sentence  in 
Step  3  and  to  remove  Step  4.  Finally,  additional  trial 
insertions  could  be  attempted  beyond  simply  replacing 
a  reactant  molecule  with  a  product  molecule.  This  may 
improve  the  acceptance  rate  in  cases  where  the 
molecular  model  character  of  the  corresponding  product 
and  reactant  are  markedly  dissimilar.  If  additional  trial 
insertions  are  attempted,  then  Ay  ^  1  in  Steps  3  and  4. 

2.5.  Choice  of  A,- 

The  efficiency  of  the  CB-RxMC  algorithm  depends  on 
the  choice  of  kh  the  number  of  trial  insertions  for  a 
given  product  molecule.  In  principle,  k,  can  be  chosen 
without  restrictions;  however,  there  is  a  computational 
expense  to  be  paid  for  increasing  the  number  of 
attempted  insertions.  The  optimal  choice  of  A,-  will 
allow  for  sampling  of  a  wide  portion  of  phase  space  but 
in  a  cost  effective  manner.  A  straightforward  approach 
for  selecting  A,-  is  to  adjust  its  value  based  on  a  desired 
acceptance  ratio  for  each  reaction.  A,-  could  then  be 
adjusted  automatically  during  the  simulation,  analogous 
to  the  adjustments  which  are  typically  made  to  the 
maximum  allowable  changes  for  particle  displacements 
and  volume  changes  [25,  26].  For  some  systems, 
however,  as  k,  is  increased,  eventually  the  computational 
cost  of  additional  trial  insertions  will  outweigh  the 
benefits  of  an  increased  acceptance  ratio.  In  the  fol¬ 
lowing,  we  propose  some  simple  guidelines  for  optimiz¬ 
ing  the  CB-RxMC  algorithm  with  respect  to  A,-.  The 
approach  is  based  on  previous  work  that  optimized  a 
configurational-bias  Monte  Carlo  method  applied  to 
chain  molecules  [27]. 

We  approximate  the  efficiency  of  the  CB-RxMC 
scheme  to  be:  (a)  proportional  to  the  probability  that  a 
given  trial  configuration  is  successfully  generated;  and 
(b)  inversely  proportional  to  the  computational  cost  of 
generating  a  configuration.  We  can  write  then  that 

( acceptance ) 

Efficiency  —  - — - - — -  (18) 

{cost) 

where  ( acceptance )  is  the  average  acceptance  ratio 
defined  as  the  number  of  accepted  reaction  steps  over 
the  number  of  attempted  reaction  steps  for  the  entire 
simulation  run,  while  {cost)  is  the  average  computa¬ 
tional  cost  of  generating  an  attempted  reaction.  In 
practice,  a  series  of  short  simulations  can  be  performed 
for  various  A/s  to  determine  the  Efficiency  as  a  function 
of  A;.  The  maximum  in  the  Efficiency  vs.  A,-  plot  will 
then  provide  the  optimal  value  for  A,-.  As  the  number 


of  trial  insertions  increases,  the  value  of  {cost)  rises. 
From  equation  (18),  we  expect  that  for  systems  where 
the  probability  of  successful  insertions  will  be  high,  e.g. 
low  density  systems,  the  optimal  value  of  A,-  will  increase 
monotonically  with  increasing  A,.  However  for  high 
density  systems  in  which  {acceptance)  marginally 
increases  with  A,-,  the  optimal  value  of  k,  will  shift  to 
lower  values. 

Note  that  for  multiple  reaction  systems  the  value  of 
A,  is  not  required  to  be  the  same  for  all  reactions,  but 
it  does  need  to  be  the  same  for  each  reactant-to-product 
replacement  pair  for  a  particular  reaction.  Also,  the 
choice  of  A,-  will  strongly  depend  on  the  state  point 
considered  as  well  as  the  details  of  the  molecular  models. 
Typical  values  of  A,-  range  from  2-200  attempted 
insertions  since  reaction  acceptance  ratios  can  vary 
widely. 

3.  Application 

We  illustrate  the  method  using  the  reacting  system, 
2H20  +  C2H4  o  2CO  +  4H2.  The  forward  reaction 
requires  the  insertion  of  three  product  molecules  into 
the  simulation  box  and  therefore  provides  a  stringent 
test  for  the  CB-RxMC  algorithm. 

3.1.  Molecular  models  and  computational  details 

The  molecular  species  C2H4,  CO,  and  H2  are  modelled 
as  single  spherical  particles  interacting  through  the 
exponential-6  potential  where  electrostatic  contributions 
are  ignored  [28].  The  water  model  consists  of  an 
exponential-6  group  at  the  oxygen  centre  and  a  set  of 
three  fixed-point  charges  [29], 

Unlike  interactions  between  species  are  approximated 
by  the  Lorentz-Berthelot  mixing  rules  [30].  A  spherical 
cutoff  of  1.05  nm  was  applied  for  the  particle-particle 
interactions  with  standard  corrections  for  this  trunca¬ 
tion  added  [26].  Long-range  electrostatic  interactions 
were  determined  by  using  the  Ewald  summation 
method  [26].  Applying  a  prescription  for  optimizing 
the  Ewald  parameters  [31,  32]  led  to  approximately 
1000-1200  wave  vectors  used  in  the  k-space  sum  where 
the  same  truncation  distance  as  above  was  applied  for 
the  real  space  sum.  The  ideal-gas  partition  functions  for 
H20,  C2H4,  CO,  and  H2  that  are  required  in  the  RxMC 
method  were  taken  from  standard  thermochemical 
reference  data  [22]. 

Simulations  were  performed  in  steps,  where  a  step 
(chosen  with  equal  probability)  was  either  a  molecule 
displacement  and  rotation  (for  H20  molecules  only), 
a  forward  reaction  step,  or  a  reverse  reaction  step. 
A  change  in  the  simulation  cell  volume  was  attempted 
every  3000  steps.  The  maximum  allowable  changes  for 
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the  displacement,  rotation,  and  volume  change  steps 
were  all  adjusted  to  achieve  an  acceptance  fraction 
of  0.4.  Simulations  were  equilibrated  for  2-3  x  106  steps 
after  which  averages  of  the  desired  quantities  were 
taken  over  5  x  106  steps.  Reported  uncertainties  were 
estimated  using  the  method  of  block  averages  [25]. 

3.2.  Results 

A  series  of  constant-pressure  RxMC  simulations  for  the 
2H20  +  C2H4  2CO  +  4H2  reaction  were  performed 
at  a  temperature  of  900  K.  The  initial  configuration  for 
each  simulation  was  250  H20  molecules  and  350  C2H4 
molecules  placed  randomly  on  a  face-centred-cubic 
lattice.  Several  pressures  were  considered  (100,  300  and 
500  MPa)  that  span  a  range  of  densities.  Further,  a  wide 
range  of  kj  values  were  considered:  1,  5,  10,  20,  40, 
80,  and  120.  Simulations  where  kj  —  1,  represent  conven¬ 
tional  sampling  of  reactions  steps  in  the  RxMC  method. 

CB-RxMC  simulation  results  are  presented  in 
table  1.  The  average  configurational  energy,  U,  and 
specific  volume,  v,  are  reported  along  with  the  average 


number  of  molecules  for  each  species,  A(.  As  expected, 
for  all  values  of  kt  the  calculated  quantities  are  within 
estimated  uncertainties.  Effects  of  increasing  kj  can  be 
seen  by  considering  the  changes  in  the  acceptance  ratio 
which  is  also  reported  in  table  1.  The  definition  of  the 
acceptance  ratio  is  the  same  as  that  used  in  equation 
(18).  (Note  that  at  equilibrium  the  forward  and  reverse 
reaction  acceptance  ratios  will  be  equal,  therefore  only 
one  value  is  reported.)  Table  1  shows  that  as  the  number 
of  trial  insertions  increases,  the  acceptance  ratio  also 
increases.  The  acceptance  ratios  for  the  three  pressures 
considered  improve  by  factors  of  about  3,  1 1  and  23, 
respectively. 

Next,  we  can  determine  an  estimate  for  the  optimal 
value  of  k,  using  the  guidelines  presented  in  section  2.5. 
In  order  to  dimensionalize  {cost)  in  equation  (18),  we 
have  introduced  a  unit  of  computational  cost  to  be  the 
time  needed  to  compute  equation  (16)  when  kj—  1  for 
all  species.  A  plot  of  the  Efficiency  as  a  function  of  kj  is 
given  in  figure  1 ,  where  the  plotted  values  are  also  given 
in  table  2.  Ensemble  averages  of  ( acceptance )  and  {cost) 
were  determined  from  the  entire  simulation  run.  All 


Table  1.  CB-RxMC  simulation  results  for  the  2H20  +  C2H4 -o  2CO  ±  4H2  reaction.3 

N,b 

-  Acceptance 


kj 

U  [kJ/mol] 

v  [cm3/g] 

h2o 

CO 

h2 

c2h4 

ratio 

P  = 

1 

100  MPa 

-2.930  (0.031) 

4.428  (0.011) 

201.3  (0.3) 

48.7  (0.3) 

96.5  (0.6) 

325.7  (0.2) 

0.1329 

5 

-2.929  (0.032) 

4.422  (0.017) 

201.2  (0.2) 

48.7  (0.2) 

97.5  (0.5) 

325.6  (0.1) 

0.3445 

10 

-2.937  (0.023) 

4.427  (0.015) 

201.3  (0.1) 

48.7  (0.2) 

97.3  (0.3) 

325.7  (0.1) 

0.3718 

20 

-2.930  (0.035) 

4.426  (0.014) 

201.3  (0.3) 

48.7  (0.3) 

97.4  (0.5) 

325.6  (0.1) 

0.3870 

40 

-2.950  (0.033) 

4.422  (0.018) 

201.4  (0.3) 

48.6  (0.3) 

97.1  (0.5) 

325.7  (0.1) 

0.3939 

80 

-2.933  (0.039) 

4.423  (0.020) 

201.3  (0.3) 

48.7  (0.3) 

97.4  (0.6) 

325.7  (0.2) 

0.3989 

120 

-2.922  (0.031) 

4.429  (0.016) 

201.3  (0.2) 

48.7  (0.2) 

97.4  (0.5) 

325.6  (0.1) 

0.4017 

P  = 

1 

300  MPa 

-6.526  (0.128) 

2.169  (0.006) 

232.2  (0.3) 

17.7  (0.3) 

34.5  (0.6) 

341.2  (0.2) 

0.0141 

5 

-6.713  (0.147) 

2.157  (0.012) 

232.6  (0.3) 

17.4  (0.3) 

34.7  (0.5) 

341.3  (0.1) 

0.0844 

10 

-6.494  (0.117) 

2.170  (0.008) 

232.3  (0.3) 

17.7  (0.3) 

35.5  (0.5) 

341.1  (0.1) 

0.1212 

20 

-6.553  (0.127) 

2.167  (0.008) 

232.3  (0.3) 

17.7  (0.3) 

35.3  (0.6) 

341.2  (0.2) 

0.1361 

40 

-6.518  (0.110) 

2.169  (0.007) 

232.3  (0.3) 

17.7  (0.3) 

35.4  (0.6) 

341.1  (0.2) 

0.1461 

80 

-6.487  (0.084) 

2.170  (0.005) 

232.2  (0.2) 

17.8  (0.2) 

35.5  (0.4) 

341.1  (0.1) 

0.1508 

120 

-6.368  (0.130) 

2.176  (0.008) 

231.9  (0.4) 

18.1  (0.3) 

35.1  (0.5) 

341.0  (0.1) 

0.1565 

P  = 

1 

500  MPa 

-7.819  (0.089) 

1.755  (0.005) 

239.5  (0.3) 

10.5  (0.3) 

20.9  (0.5) 

344.7  (0.1) 

0.0037 

5 

-8.024  (0.242) 

1.746  (0.006) 

239.9  (0.2) 

10.2  (0.2) 

20.4  (0.5) 

344.9  (0.1) 

0.0321 

10 

-7.818  (0.142) 

1.753  (0.003) 

239.5  (0.2) 

10.4  (0.2) 

20.9  (0.4) 

344.8  (0.1) 

0.0549 

20 

-7.758  (0.080) 

1.756  (0.003) 

239.5  (0.2) 

10.5  (0.2) 

21.0  (0.3) 

344.7  (0.1) 

0.0707 

40 

-7.742  (0.130) 

1.757  (0.005) 

239.4  (0.2) 

10.6  (0.2) 

21.1  (0.3) 

344.7  (0.1) 

0.0814 

80 

-7.834  (0.129) 

1.756  (0.005) 

239.5  (0.2) 

10.5  (0.2) 

21.0  (0.4) 

344.8  (0.1) 

0.0844 

120 

-7.859  (0.117) 

1.753  (0.004) 

239.5  (0.2) 

10.5  (0.2) 

20.9  (0.3) 

344.8  (0.1) 

0.0852 

3Uncertainty  in  units  of  the  last  decimal  digit  is  given  in  parentheses:  —2.930  (0.031)  implies  —2.930  ±0.031. 
bNumber  of  molecules  of  species  i. 
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Figure  1.  Efficiency  of  the  CB-RxMC  algorithm  as  a  function  of  the  number  of  trial  insertions  for  the 
2H20  +  C2H4  -o-  2CO  +  4H2  reaction  at  the  three  pressures  considered.  Line  is  shown  as  a  guide  to  the  eye  only. 


Table  2.  Efficiency  at  various  kt  values  for  the 
2H20  +  C2H4  2CO  +  4H2  reaction. 


ki 

Efficiency11 

100  MPa 

300  MPa 

500  MPa 

1 

0.133 

0.014 

0.0037 

5 

0.343 

0.084 

0.0319 

10 

0.364 

0.119 

0.0538 

20 

0.375 

0.132 

0.0684 

40 

0.369 

0.137 

0.0762 

80 

0.349 

0.132 

0.0739 

120 

0.332 

0.129 

0.0704 

Efficiency  determined  from  equation  (18). 


three  curves  exhibit  similar  behaviour  although  the 
behaviour  is  most  pronounced  for  the  100  MPa  case. 
Comparing  the  three  curves,  for  all  values  of  kt  the 
efficiency  decreases  as  the  pressure  (and  corresponding 
density)  of  the  system  increases.  This  is  expected  since 
the  space  available  for  product  molecule  insertion  is  also 
reduced  with  increasing  pressure.  The  maximum  for  the 
100  MPa  system  is  at  approximately  k,  —  20,  while 
the  maximum  for  the  300  and  500  MPa  systems  is  at 
approximately  =  40.  However,  significant  gains  in 
the  efficiency  of  the  CB-RxMC  algorithm  are  found 
for  rather  small  values  of  kt  just  beyond  conventional 
sampling  (k,  =  1),  i.e.  k,  =  5, 10.  Beyond  the  maximum 


values  of  /c,  for  each  pressure  considered,  the  efficiency 
begins  to  decrease.  This  can  more  clearly  be  seen  by 
considering  table  2.  For  the  300  and  500  MPa  cases  it 
is  evident  that  attempting  more  than  40  trial  insertions 
per  reaction  step  is  not  beneficial.  Moreover,  given  the 
flatness  of  the  curves  at  high  k,  as  well  as  the  steep  rise  at 
low  kb  the  most  practical  value  for  k,  is  approximately 
15-20  trial  insertions  per  attempted  reaction.  For  the 
100  MPa  case,  the  optimal  value  of  kf  is  more  obvious, 
i.e.  kf—  20. 


4.  Conclusions 

A  cavity-bias  algorithm  was  presented  that  increases 
the  acceptance  rate  of  reaction  steps  in  the  reaction 
ensemble  Monte  Carlo  method.  For  some  cases,  the 
CB-RxMC  method  is  shown  to  increase  acceptance 
rates  for  the  reaction  step  by  a  factor  greater  than  20. 
CB-RxMC  enhances  phase  space  sampling  but  more 
critically  allows  for  the  simulation  of  systems  not 
feasible  with  conventional  random  sampling.  CB-RxMC 
is  an  extension  of  other  cavity-bias  sampling  techniques 
used  in  the  grand  canonical  ensemble  MC  method  [18]. 
The  main  bottleneck  of  an  RxMC  simulation  is 
the  insertion  of  product  molecules  into  the  simulation 
box  during  a  reaction  step.  In  a  biased  fashion,  the 
CB-RxMC  method  generates  energetically  favourable 
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positions  for  these  molecules  by  searching  for  cavities  in 
the  system.  The  bias  introduced  when  generating 
positions  is  then  corrected  for  in  the  acceptance  criteria. 
For  a  given  state  point,  the  efficiency  of  the  CB-RxMC 
technique  depends  on  the  number  of  trial  insertions  (kj) 
for  each  attempted  reaction  step.  A  tradeoff  exists 
between  the  additional  computational  expense  of  gen¬ 
erating  kj  configurations  for  a  single  attempted  reaction 
and  overall  phase  space  sampling.  A  simple  approach 
was  suggested  for  optimizing  the  value  of  kj. 

With  minor  modifications,  the  algorithm  presented 
here  could  include  orientational-bias  techniques  for 
polar  molecules  as  well  as  biased  sampling  techniques 
for  inserting  long-chain  molecules  [25  and  references 
therein]. 
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